Ce tutoriel montre quelques exemples de traitements géographiques à partir de la base Aspe, afin de :
Le chargement du package {aspe} est supposées déjà réalisé (voir ce support), de même que l’importation initiale des données.
library(aspe)
library(tidyverse)
library(mapview)
load(file = "processed_data/toutes_tables_aspe_sauf_mei.RData")
La table ref_type_projection fait correspondre un code (typ_id) à différentes manières de désigner le système de projection (le CRS). Ici nous utiliserons en particulier le libellé (typ_libelle_sandre) et le code EPSG de chacun des CRS.
Pour des traitements géographiques, il est nécessaire que tous les objets manipulés soient dans un même CRS.
Quelle est la fréquence de chaque CRS parmi les stations et les points de prélèvement ?
crs_stations <- station %>%
rename(typ_id = sta_typ_id) %>%
left_join(y = ref_type_projection) %>%
pull(typ_libelle_sandre) %>%
as.character() %>%
table() %>%
as.data.frame() %>%
purrr::set_names(c("CRS", "Nombre de stations"))
crs_points <- point_prelevement %>%
rename(typ_id = pop_typ_id) %>%
left_join(y = ref_type_projection) %>%
pull(typ_libelle_sandre) %>%
as.character() %>%
table() %>%
as.data.frame() %>%
purrr::set_names(c("CRS", "Nombre de points"))
crs_tot <- full_join(x = crs_points, y = crs_stations)
DT::datatable(crs_tot, width = 600, rownames = FALSE)
Concernant les CRS utilisés pour les stations et les points de prélèvement, la table ref_type_projection contient les données suivantes :
ref_type_projection %>%
filter(typ_libelle_sandre %in% (crs_tot %>% pull(CRS))) %>%
as.data.frame() %>%
DT::datatable()
Il apparaît que pour le Lambert II étendu, le code EPSG (27572) est manquant. Il faut donc compléter la table ref_type_projection.
ref_type_projection <- ref_type_projection %>%
mutate(typ_code_epsg = ifelse((is.na(typ_code_epsg) & typ_libelle_sandre == "Lambert II Etendu"),
yes = 27572,
no = typ_code_epsg))
Vérification :
ref_type_projection %>%
filter(typ_libelle_sandre %in% (crs_tot %>% pull(CRS))) %>%
as.data.frame() %>%
DT::datatable()
Comme la base comprend les outre-mers, il faut utiliser un CRS mondial, ce qui n’est pas le cas des systèmes Lambert. On choisit donc le WGS84 qui est utilisé par les GPS, OpenStreetMaps, etc.
station des codes EPSG de chacune des stationsstation <- station %>%
left_join(y = ref_type_projection,
by = c("sta_typ_id" = "typ_id"))
ATTENTION, cette opération prend une dizaine de minutes environ.
coords_wgs84 <- geo_convertir_coords_df(df = station,
var_x = "sta_coordonnees_x",
var_y = "sta_coordonnees_y",
var_crs_initial = "typ_code_epsg",
crs_sortie = 4326) %>%
rename(x_wgs84 = X, y_wgs84 = Y)
Comme la reprojection est une opération chronophage, il est intéressant d’en sauvegarder le résultat.
# save(coords_wgs84, file = "processed_data/sta_coords_wgs84.RData")
Ajout des coordonnées WGS84 et suppression des colonnes qui ne serviront plus :
station <- station %>%
bind_cols(coords_wgs84) %>%
select(-(sta_geometrie:typ_code_epsg))
Vérification que les colonnes correspondant aux coordonnées en WGS84 ont bien été créées.
names(station)
## [1] "sta_id" "sta_code_sandre"
## [3] "sta_libelle_sandre" "sta_enh_id"
## [5] "sta_eligibilite_ipr_iprplus" "sta_com_code_insee"
## [7] "sta_point_km_aval" "sta_localisation_precise"
## [9] "sta_code_national_masse_eau" "x_wgs84"
## [11] "y_wgs84"
A ce stade, notre objet station est un simple tableau de données (dataframe) avec des colonnes correspondant aux coordonnées. Pour effectuer une sélection spatiale, il faut manipuler des objets spatiaux comme en SIG. Ces manipulations sont relativement simples grâce au package sf.
Il s’agit de créer l’homologue d’une couche SIG des stations.
station_geo <- station %>%
sf::st_as_sf(coords = c("x_wgs84", "y_wgs84"), crs = 4326)
Visualisation.
mapview::mapview(station_geo)